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^ ■ Abstract 
> 

Q ' Time-symmetric integration schemes share with symplectic schemes the property 

■ that their energy errors show a much better behavior than is the case for generic 
integration schemes. Allowing adaptive time steps typically leads to a loss of sym- 
plecticity. In contrast, time symmetry can be easily maintained, at least for a con- 
tinuous choice of time step size. In large-scale N-body simulations, however, one 

, often uses block time steps, where all time steps are forced to take on values as 

I powers of two. This greatly facilitates parallelization, and hence code efficiency. 

■ Straightforward implementation of time-symmetry, translated to block time steps, 
faces significant hurdles. For example, iteration can lead to oscillatory behavior, and 
even when such behavior is suppressed, energy errors show a linear drift in time. 
We present an approach that circumvents these problems. 
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tj ■ 1 Introduction 



For long-term N-body simulations, it is essential that the drift in the values 
of conserved quantities is kept to a minimum. The total energy is often used 
as an indicator of such a drift. During the last fifteen years, two approaches 
have been put forward to improve numerical conservation of energy and other 
theoretically conserved quantities: symplectic integration schemes, where the 
simulated system is guaranteed to follow a slightly perturbed Hamiltonian 
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system, and time-symmetric integration schemes, where the simulated system 
follows the same trajectory in phase space, when run backwards or forwards 



In both cases, for symplectic as well as for time-symmetric schemes, the intro- 
duction of adaptive time steps tends to destroy the desired properties. Sym- 
plectic schemes are perturbed to different Hamiltonians at different choices of 
time step length, and therefore lose their global symplecticity. Time-symmetric 
schemes typically determine their time step length at the beginning of a step, 
which implies tha t running a step backwards gives a slightly different length 
for that step. See iLeimkuhler fc Reich (Hooi) for a recent review of various 



attempts to remedy this situation. 

In practice, many large-scale simulations in stellar dynamics use a block time- 
step approach, where the on ly allowed values for the time step length are 
powers of two (Aarseth, 2003[ ). The name derives from the fact that, with this 



recipe, many particles will share the same step size, which implies that their 
orbit integration can be performed in parallel. An added benefit, in the case 
of individual time steps, is the fact that block time steps allow one to predict 
the positions of all particles only once per block time step, rather than sepa- 
rately for each particle that needs to be moved forward. Since parallelization is 
rapidly becoming essential for any major simulation, we explore in this paper 
the possibility to extend time symmetry to the use of block time steps. 

In Section 2, we analyze some of the problems that occur when applying 
existing methods to the case of block time steps, and we offer a novel solution, 
with a truly time symmetric choice of time step, with the restriction that we 
only allow changes of a factor two in the direction of increasing and decreasing 
the time step. In Section 3, we present numerical tests of our new scheme, 
which show the superiority of our approach over various alternatives. Section 
4 sums up. 



2 Block Time Steps 



2. 1 Implicit Iterative Time Symmetrization 



It is surprisingly easy to introduce a time-symmetric version for any adaptive 
self-starting integration scheme. Let ^ = (r, v) be the 2A^-dimensional phase 
space vector for a system with degrees of freedom, and let /(^j, 5ti) be the 
operator that maps the phase space vector of the system at time ti to a new 
phase space vector at time tj+i = ti+5ti. Any choice of self-starting integration 
scheme, together with a recipe to determine the next time step 6ti, at time tj 
and phase space value C,iiti), defines the precise form of /(^j, Sti). 
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The re cipe for making any such scheme time-symmetric was given bv lHut. Making fc McMillan 
fjl99.5h . as: 



6+1 = f{^i,Sti), 



6t 



where /i(^) can be any time step criterion. Note that this recipe leads to an 
implicit integration scheme, which can be solved most easily through iteration. 
In practice, one or two iterations suffice to get excellent accuracy, but at 
the cost of doubling or tripling the number of force calculations that need 
to be perfor nied. Extensions of t his i mplicit symmetri zation idea have been 
presented bv lFunato et al. (|l99fit) and lHut et all (|l997l) . 



Since we will need to inspect the idea of iteration below in more detail, let us 
write out the process here explicitly. We start with the given state and the 
implicit equation for ^j+i of the form 



e.+i = /(6,5t.(6,6+i)) (2) 

The ffist guess for ^j+i is 

d?i = /(6,'5t.(6,6)) (3) 

and we can consider this as our zeroth-order iteration. With this guess in hand, 
we can now start to iterate, finding 



d+\ = /(6,5t.(6,d?i)) (4) 

as our ffist-order iteration. This will already be much closer to the final value, 
as long as the time steps are small enough and the function 6ti does not 
fluctuate too rapidly. In general, the /c*^ iteration will yield a value for ^j+i of 



&\ = fi^^,StM^,^tl'^)) (5) 

We will now consider the application of these techniques to block time steps. 
For the purpose of illustrating the use of block time steps, it will suffice to 
use the leapfrog scheme (also known as the Verlet-Stormer-Delambre scheme, 
according to the authors who rediscovered this scheme at roughly century-long 
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intervals), which we present here in a self-starting, but still time-symmetric 
form: 



Tj+i = + Vi6t + ai{6tf/2 
Vi+i = Vi + (tti + ai+i)5t/2 



(6) 



2.2 Flip-Flop Problem 



To start with, we apply the recipe of lHut. Makino fc McMilk^ (|l99,^ to block 
time steps. Let us define a block time step at level n as having a length: 



Atr 



Atl 

On-l ■ 



(7) 



where Ati is the maximum time step length. Starting with the continuum 
choice of 



Sir 



+1/ 



(8) 



we now force each time step to take on the block value 5ti = At„ for the 
smallest n value that obeys the condition Atn < Stc^i. In more formal terms, 
6ti = 6ti{6tc,i) = Atn for the unique n value for which 



n = min < k 

k>l 



1 ^ hj^i) + /ife+i) 



2k-i 



(9) 



The problem with this approach is that we are no longer guaranteed to find 
convergence for our iteration process, as can be seen from the following exam- 
ple. Let h{C,i) = 0.502 and let the time derivative of h{C,i{t)) along the orbit 
be {d/dt)h{^i(t)) = — O.OL We then get the following results for our attempt 
at iteration. 



= /(6,5t.(0.502)) = /(e,,0.5) 



(10) 
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Sti{[0.502 + (0.502 + 0.5 * (-0.01))]/2)) 
(5ti([0.502 + 0.497]/2)) 



,11 



(11) 




/(e 
/(e 



5^.(6, e^+i)) = /(6,5t.([/^(e.) + Meli)]/2)) 

(5ti([0.502 + (0.502 + 0.25 * (-0.01))]/2)) 
5ti([0.502 + 0.4995]/2)) 
5ti(0.50075)) = /(ei,0.5) 



i5 



.1) 



(12) 



And from here on, = /(^j, 0.25) for every odd value of k and ^l_Ji = 
f{^i, 0.5) for every even value of k: the process of iteration will never converge. 

Under realistic conditions, for slowly varying h functions and small time steps, 
this flip-flop behavior will not occur often, but it will occur sometimes, for a 
non-negligible fraction of the time. We can see this already from the above 
example: for a linear decrease in the h function of [d / dt)h{^i{t)) = —0.01, we 
will get flip-flopping not only for h{^i) — 0.502 but for any value in the finite 
range 0.50125 < h{^i) < 0.5025. 

Since iteration converges correctly over the rest of the interval 0.5 < h{^i) < 1, 
we conclude that in this particular case flip-flopping occurs about one quarter 
of one percent of the time, over this interval. This is far too frequent to be 
negligible in a realistic situation. 

Clearly, a straightforward extension of the implicit iterative time symmetriza- 
tion approach does not work for block time steps, because iteration does not 
converge. We have to add some feature, in some way. Our flrst attempt at a 
solution is to take the smallest of the two values in a flip-flop situation. 

2.3 Flip-Flop Resolution 

The most straightforward solution of the flip-flop dilemma is like cutting the 
Gordian knot: we just take the lowest value of the two alternate states. The 
drawback of this solution is that in general we need at least two iterations 

for each time step, to make sure that we have spotted, and then correctly 
treated, all flip-flop situation. In general, it is only at the third iteration that it 
becomes obvious that a flip-flop is occurring. To see this, consider the previous 
example with a starting value of h{^i) ~ 0.501. In that case we will get ^f^^ — 



/({i,0.5) and 4>i = /(Cj)0-25), just as when we started with h{$,i) — 0.502. 
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The difference shows up only at the second iteration, where we now find = 
/(^i, 0.25), a value that will hold for all higher iterations as well. 

The original iterative approach to time symmetry in practice already gives 
good results when we use only one iteration. This implies a penalty, in terms 
of force calculations per time steps, of a factor two compared to non-time- 
symmetric explicit integration. Now the use of flip-flop resolution will force us 
to always take at least two iterations per step, raising the penalty to become 
at least a factor of three. 

However, there is a more serious problem: there is still no guarantee that 
taking the lowest value in a flip-flop situation leads to a time-symmetric recipe. 
In fact, what is even more important, we have not yet checked whether our 
symmetric block time-step scheme is really time symmetric, in the absence of 
flip-fiop complications. 

In order to investigate these questions, let us return to the example we used 
above, but instead of a linear time derivative, let us now use a quadratic time 
derivative for the h function that gives the estimate for the time step size. 
Rather than writing a formal definition, let us just state the values, while 
shifting the time scale so that t — coincides with the particle position being 
6: 



/i(0.00) = 0.502 
/i(0.25) = 0.499 

/i(0.50) = 0.499 (13) 
When we start at time t — 0, and we integrate forwards, we find: 

d+i=/te,5t.(6,6)) = f{^^,5u{h{m 

= /(e„5t.(0.502)) = /(e„0.5) (14) 



= /(ei,5t,(0.5005)) = /(ei,0.5) (15) 

and so on: all further kth iterations will result in ^^^j^. = /(^j, 0.5). There is no 
fiip-fiop situation, when moving forward in time. 

However, when we now turn the clock backward, after taking this step of 
half a time unit, we start with the value /i(0.50) = 0.499, which leads to a 
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first step back of 5t = 0.25. Tlie end point of the first step back is t = 0.25 
witli h{0.25) = 0.499. Tlierefore, also liere tliere is no flip-flop situation: all 
iterations, while going backward, result in a time step size of 5t = 0.25. 

We have thus constructed a counter example, where forward integration would 
proceed with time step 6t = 0.5 and subsequent backward integration would 
proceed with time step dt = 0.25. Clearly, our scheme is not yet time sym- 
metric, even in the absence of a flip-flop case. 

2.4 A First Attempt at a Solution 

Let us rethink the whole procedure. The basic problem has been that the 
very first step in any of our algorithms proposed so far has not been time 
symmetric. The very first step moves forward, and leads to a newly evolved 
system at the end of the first step. Only after making such a trial integration, 
do we look back, and try to restore symmetry. However, as we have seen, the 
danger is large that this trial integration is not exhaustive: it may already go 
too far, or not far enough, and thereby it may simply overlook a type of move 
that the same algorithm would make if we would start out in the time-reversed 
direction. 

Formulating the problem in this way, immediately STiggests a solution. At any 
point in time, let us first try to make the largest step that is allowed. If that 
step turns out to be too large for our algorithm, we try a step that is half that 
size, if that step is too large still, we again half the size, and so on, until we 
find a step size that agrees with our algorithm, when evaluated in both time 
directions. A similar treatment has been described by Quin et al (1997). 

This type of approach is clearly more symmetric than what we have attempted 
so far. Instead of using information of the physical system at the starting point 
of the next integration step, we only use a mathematical criterion to find the 
largest time step size allowed at that point, and we then apply the physical 
criteria symmetrically in both directions. 

Let us give an example. If the largest time step size is chosen to be unity, then 
at time t = we start by considering this time step. Wc try, in this order 
6t = 1, 6t = 0.5, 6t = 0.25, and so on, until we find a time step for which 
integration starting in the forward direction, and integration starting in the 
backward direction, both result in the new time step being acceptable. Let us 
say that this is the case for 5t = 0.125. 

After taking this step, we are at time t = 0.125. The largest time step allowed 
at that point, forward or backward, is 6t = 0.125. Any larger time step would 
result in non-alignment of the block time steps: in the backward direction it 
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would jump over t = 0. So at this point we start by considering once more 
5t = 0.125. If that time step is too large, we try half that time step, halving 
it successively until we find a satisfactory time step size. 

Imagine that the second time step size is also St = 0.125. In that case, we land 
at t = 0.25. From there on, the maximum allowed time step size is 5t — 0.25, 
so the first try should be that size. 

In principle, this approach seems to be really time symmetric. However, there is 
a huge problem with this type of scheme, as we have just formulated it. Imagine 
the system to crawl along with time steps of, say 6t = 1/1024, and reaching 
time t = 1. Our new recipe then suggests to start by trying 6t = 1, a 1024-fold 
increase in time step! Whatever subtle physical effect it was that forced us to 
take such small time steps, is completely ignored by the mathematical recipe 
that forces us to look at such a ridiculously large time step. 

For example, in the case of stellar dynamics, a double star may force the stars 
that orbit each other to take time steps that are necessarily far shorter than 
the orbital period. Starting out with a trial step size that is far larger than 
an orbital period may or may not give spuriously safe-looking results. Clearly, 
we have to exclude such enormous jumps in time step. 

2.5 A Second Attempt at a Solution 

The simplest solution to taming sudden unphysical increases in time steps 
is to allow at most an increase of a factor two, in either the forward or the 

backward direction. This then implies that we can only allow decreases of a 
factor two, and not more than two, in either direction. The reason is that a 
decrease of a factor four in one direction in time would automatically translate 
into an increase of a factor four in the other direction. 

Note that we have to be careful with our time step criterion. If we allow time 
steps that are too large, we may encounter situations where our time step 
criterion would suggest us to shrink time steps by a factor of four, from one step 
to the other. Since our algorithm does not allow this, we can at most shrink by 
a factor of two, which may imply an unacceptably large step. However, if our 
time step criterion is sufficiently strict, allowing only reasonably small time 
steps too start with, it will be able to resolve the gradients in the criterion in 
such as way as to handle all changes gracefully through halving and doubling. 

When we apply this restriction to the scheme outlined in the previous subsec- 
tion, we arrive at the following compact algorithm. 

First a matter of notation. Any block time step, of size St — 1/2^, connects 
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two points in time, only one of which can be written at t = Z I2'^^~^\ with Z 
an integer. Let us call that time value an even time, from the point of view of 
the given time step size, and let us call that other time value an odd time. To 
give an example, if 5t — 0.125, than t — 0, 0.25, 0.5, 0.75, 1 are all even times, 
while t = 0.125, 0.375, 0.625, 0.875 are all odd times. 

Here is our algorithm: 

When we start in a given direction in time, at a given point in time, we 
should determine the time step size of the last step made by the system. In 
that way, we can determine whether the current time is even or odd, with 
respect to that last time step. 

If the current time is odd, our one and only choice is: to continue with the 
same size time step, or to halve the time step. First, we try to continue with 
the same time step. If, upon iteration, that time step qualifies according to 
the time-symmetry criterion used before, Eq. 9, we continue to use the same 
time step size as was used in the previous time step. If not, we use half of 
the previous time step. 

If the current time step is even, we have a choice between three options 
for the new time step size: doubling the previous time step size, keeping it 
the same, or halving it. We first try the largest value, given by doubling. 
If Eq.9 shows us that this larger time step is not too large, we accept it, 
otherwise we consider keeping the time step size the same. If Eq.9 shows 
us that keeping the time step size the same is okay, we accept that choice, 
otherwise we just halve the time step, in which case no further testing is 
needed. 

Note that in this scheme, we always start with the largest possible candidate 
value for the time step size. Subsequently, we may consider smaller values, 
but the direction of consideration is always from larger to smaller, never from 
smaller to larger. This guarantees that we do not run into the flip-flop problem 
mentioned above. 



3 Numerical Tests 

We present here the results for a gravitational two-body integration. The rel- 
ative orbit of the two point masses forms an ellipse with an eccentricity of 
e = 0.99. We have chosen a time unit such that the period of the orbit is 
r= 27r. 

We have implemented four different integration schemes: 
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0) the original time-symmetric integration scheme described bv lHut. Makino fc McMillai 
(EH), where there is a continuous choice of time step size. This is the ap- 
proach described in section 2.1. We have used five iterations for each step. 



1) a block-time-step generalization, with a fixed number of iterations. This is 
the approach analyzed in section 2.2. Here, too, we chose five iterations for 
each step. 

2) a block time step generalization, with a variable number of iterations. If 
after five iterations, the fourth and the fifth iterations still give a different 
block time step size, then we choose the smallest of the two. This recipe 
avoids fiip-fiop situations. It is the approach described in section 2.3. 

The algorithm described in the next section, 2.4, we have not implemented 
here, because it is guaranteed to lead to large errors in those cases where a new 
large time step is allowed again just before pericenter passage. We therefore 
switched directly to the following section: 

3) the implementation of our favorite algorithm, where we start with a truly 
time symmetric choice of time step, with the restrictions that we only allow 
changes of a factor two in the direction of increasing and decreasing the time 
step, and that we only allow an increase of time step on the so-called even 
time boundaries. This is the approach given in section 2.5. 

In figures 1 and 2 we show the results of integrating our highly eccentric 
binary with these four integration schemes. In each case, the largest errors are 
produced by algorithm 1), smaller errors are produced by algorithm 2), and 
even smaller errors appear with algorithm 3). Finally, algorithm 0) gives the 
smallest errors. 

Figure 1 shows the energy error in the two-body integration as a function 
of time. As is generally the case for time-symmetric integration, the errors 
that occur during one orbit are far larger than the systematic error that is 
generated during a full orbit. To bring this out more clearly, figure 2 shows 
the error only one time per orbit, at apocenter, the point in the orbit where 
the two particles are separated furthest from each other, and the error is the 
smallest. 

Finally, figure 3 shows the same data as figure 2, but for a period of time that 
is ten times longer. In both figures 2 and 3, it is clear that the first two block 
time steps algorithms, 1) and 2), both show a linear drift in energy. This is 
a clear sign of the fact that they violate time symmetry. Note that in both 
figures algorithm 3) gives rise to a time dependency that looks like a random 
walk. This may well be the best that can be done with block time steps, when 
we require time symmetry. 
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Fig. 1. Energy errors for a two-body integration of a bound orbit with eccentricity 
e = 0.99. The top hnc with highest slope corresponds to algorithm 1, the line with 
intermediate slope corresponds to algorithm 2, and below those the two lines for 
algorithms and 3 are indistinguishable in this figure. 

4 Summary 



To sum up, we have succeeded in constructing an algorithm for time sym- 
metrizing block time steps that does not show a linear growth of energy errors. 
As far as we know, this is the first such algorithm that has been discovered. 
We expect this algorithm to have practical value for a wide range of large-scale 



11 



0.01 



AE 



5x10" 












2000 4000 6000 



5000 



10' 



Fig. 2. The energy errors at apocenter. The four hnes, from top to bottom, corre- 
spond to algorithms 1, 2, 3, and 0. 

parallel N-body simulations. 
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sis, available as "www-theorie.physik.unizh.ch/~stadel/downloads/thesis.ps". 
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Fig. 3. Same as Fig. 2, but for a duration that is ten times longer. 
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